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INTRODUCTION 

Forced  convection  heat  transfer  to  fluids  flowing  in  conduits  has  been 
of  interest  in  chemical  engineering  practice  for  many  years.   Heat  exchangers 
and  cooling  coils,  two  very  common  devices,  both  transfer  heat,  in  part  if 
not  totally,  by  forced  convection.   In  forced  convection  heat  is  transferred 
by  simultaneous  conduction  and  fluid  flow,  the  pattern  of  which  is  determined 
by  an  external  force  as  might  be  exerted  by  a  pump  impeller. 

The  rate  of  heat  transfer  in  conduits  is  not  uniform,  but  is  greater 
near  the  entrance  where  the  velocity  and  temperature  profiles  are  developing. 
The  nature  and  rate  of  these  developing  profiles  depend  primarily  on  the 
Reynolds  number  as  well  as  the  Prandtl  number.   A  majority  of  the  published 
annular  heat  transfer  work  assumes  uniform  temperature  and  fully  developed 
velocity  profile  at  the  entrance  (2=0)  and  therefore  results  in  reasonably 
simple  governing  equations.   In  the  present  analysis,  the  rate  of  heat  trans- 
fer is  calculated  for  laminar  flow  in  an  annulus  with  simultaneously  develop- 
ing velocity  and  temperature  profiles.   The  boundary  layer  assumptions  have 
been  made  so  that  the  work  is  confined  to  a  region  between  creeping,  or  very 
slow  flow,  and  turbulent  flow. 

It  is  well  known  from  boundary  layer  theory  that  upon  entering  a  con- 
duit a  fluid  undergoes  a  velocity  development  during  its  course  through  the 
conduit.   Longitudinally  growing  boundary  layers  are  formed  on  each  of  the 
boundary  surfaces.   The  thicknesses  of  these  layers  increase  as  the  fluid 
progresses  downstream  until  a  point  is  reached  where  the  boundary  layer  on 
one  surface  intercepts  that  from  the  opposite  surface.   From  this  point 


onward  (at  a  certain  value  of  z) ,  the  velocity  profile  will  remain  unchanged 
with  further  increases  in  the  downstream  distance  and  is  then  considered  to 
be  fully  developed. 

The  temperature  profile  develops  more  or  less  in  the  same  manner  except 
it  should  be  mentioned  that  it  depends  to  a  great  extent  on  the  velocity 
profile.   This  influence  is  indicated  by  the  presence  of  the  velocity  terms 
in  the  differential  energy  equation.   On  the  other  hand,  it  is  not  so  appar- 
ent that  the  velocity  profile  is  influenced  by  the  temperature  profile.   In 
fact,  the  momentum  equation  does  not  contain  temperature  explicitly,  but  it 
does  contain  temperature-dependent  terms  such  as  viscosity. 

The  determination  of  the  velocity  and  temperature  profiles  involves 
solution  of  the  differential  mass,  energy  and  momentum  balances.   As  is 
often  the  case,  analysis  of  systems  of  engineering  interest  results  in  very 
complex  partial  differential  equations,  analytical  solutions  of  which  exist 
only  for  simplified  cases.   In  ether  words,  assumptions  have  to  be  made  to 
reduce  the  problem  to  a  tractable  level.   In  this  study,  a  numerical  approach 
has  been  used  exclusively.   Careful  assumptions  have  been  made  to  render  the 
problem  soluble,  but  without  oversimplication. 

It  is  assumed  herein  that  the  fluid  properties  are  constant  and  that 
the  velocity  and  temperature  are  uniform  at  the  entrance  cross  section.   The 
study  is  also  restricted  to  the  consideration  of  incompressible,  laminar, 
Newtonian  flow  with  negligible  axial  conduction  and  viscous  dissipation. 
Two  distinct  problems  with  various  values  of  the  Prandtl  number  and  of  the 
ratio  of  che  inner  and  outer  radii  are  considered: 

(i)   For  z  >  0,  constant  wall  temperature  at  the  inner  wall  and  outer 
wall  insulated 


(ii)   For  z  >  0,  uniform  heat  flux  at  the  inner  wall  and  insulation 
at  the  outer  wall. 

The  assumption  of  constant  physical  properties  avoids  coupling  of  the  momen- 
tum equation  to  the  energy  equation  and  therefore  permits  determination  of 
the  velocity  profile  independently.   Once  the  flow  pattern  is  known,  the 
temperature  distribution  can  be  calculated.  ^ 

It  is  the  purpose  of  this  work  to  study  the  variation  of  the  Nusselt 
number  with  the  axial  distance  from  the  inlet  (z=0) .   Unfortunately,  due  to 
the  non-linearity  of  the  equation  of  motion,  there  is  no  known  analytical 
method  of  solution.   Two  approximate  methods  involving  linearizing  the  mo- 
mentum equation  have  been  employed  by  previous  investigators  of  this  problem. 
The  fluid  is  first  assumed  to  enter  the  conduit  with  a  constant  velocity 
parallel  to  its  axis.   With  this  initial  velocity,  the  momentum  equation  is 
linearized  either  by  Targ's  approximation  or  by  applying  the  technique  of 
Langhaar.   The  resulting  simplified  equation  can  then  be  solved,  often  by 
analytical  methods. 

No  linearization  has  been  attempted  in  this  study.   Instead,  the  finite 
difference  technique  has  been  adopted  and  the  non-linear  terms  in  the  equa- 
tion of  motion  retained.   In  this  method,  the  partial  derivatives  are  re- 
placed by  difference,  quotients  in  the  independent  variables  and  the  result 
used  as  an  approximation  of  the  derivatives.   The  partial  differential  equa- 
tions are  then  reduced  to  sets  of  simultaneous  algebraic  equations  which  can 
be  solved  with  the  aid  of  high-speed  digital  computers. 

As  a  consequence  of  the  boundary  layer  assumptions  the  Reynolds  number 
no  longer  appears  as  a  parameter.   It  should  be  borne  in  mind,  however, 
that  the  assumptions  made  in  obtaining  the  boundary  layer  equations  for 


this  system  are  satisfied  with  an  increasing  accuracy  as  the  Reynolds  number 
increases.   Unlike  a  study  with  fully  developed  velocity  profile,  the  Prandtl 
number  becomes  of  great  significance  when  both  profiles  are  developing  simul- 
taneously.  Being  the  ratio  of  momentum  and  thermal  diffusivity,  it  provides 
a  measure  of  the  relative  rate  of  formation  of  the  momentum  and  thermal 
boundary  layers.   At  a  Prandtl  number  of  10,  for  example,  the  rate  of  momen- 
tum diffusion  is  greater  than  that  of  thermal  diffusion,  with  the  result 
that  the  velocity  profile  approaches  its  fully-developed  pattern  more  rapidly 
than  does  the  temperature  profile.   On  the  contrary,  for  a  small  Prandtl 
number,  say  0.01,  the  temperature  profile  is  established  much  faster.   For 
Prandtl  numbers  near  unity,  both  velocity  and  temperature  develop  at  a  simi- 
lar rate.   In  short,  the  solutions  for  uniform  velocity  profile,  Pr  =  0,  and 
for  fully-developed  velocity  profile,  Pr  =  °°,  correspond  to  the  upper  and 
lower  limits. 

Another  important  parameter  is  the  ratio  of  the  inner  to  the  outer 
radius  of  the  annulus.   A  ratio  of  zero  gives  tubular  flow  while  a  ratio  of 
unity  corresponds  to  the  parallel-plate  flow.   These  two  cases,  which  repre- 
sent the  upper  and  lower  limits  of  annular  flow,  have  long  been  studied  by 
many  investigators.   In  fact,  if  the  limit  of  the  annular  flow  solution  is 
taken  as  this  parameter  tends  to  zero,  it  is  obvious  that  the  solution  will 
not  approach  that  of  tubular  flow,  because  for  annular  flow,  an  entirely 
different  boundary  condition  applies  at  the  inner  wall. 

In  this  analysis,  cases  for  which  the  Prandtl  number  takes  on  values 
of  0.01,  1  and  10  and  for  which  the  radius  ratio  takes  on  values  of  2.0  and 
10.0  are  studied  for  the  two  sets  of  boundary  conditions. 


CHAPTER  II 


LITERATURE  SURVEY 


The  important  problem  of  incompressible  laminar  flow  in  the  entrance 
region  of  a  circular  tube  has  been  studied  by  many  investigators  from  the 
points  of  view  of  both  momentum  and  energy  transport.   A  number  of  these 
papers  and  those  for  other  geometries  as  well  are  given  in  Ref.  (22).   In 
recent  years,  the  annulus  problem  has  also  gained  increased  significance  in 
a  number  of  engineering  developments  such  as  the  cooling  of  high-temperature 
atomic  reactors.   However,  due  to  the  complexity  of  the  governing  differen- 
tial equations,  the  problem  of  simultaneously  developing  velocity  and  tem- 
perature distributions  in  the  entrance  region  of  an  annulus  has  received 
much  less  attention.   Only  a  few  papers  pertinent  to  the  work  have  been 
published  to  date.   In  the  discussion  below,  the  literature  surveyed  will 
be  divided  into  three  categories,  namely  (i)  the  velocity  profile,  (ii)  the 
temperature  profile  and  (iii)  simultaneous  development  of  velocity  and  tem- 
perature profiles.   The  mathematical  techniques  involved  will  be  outlined 
when  necessary. 

Various  approximate  solution  methods  have  been  devised  in  the  past  for 
determining  the  laminar  flow  development  of  a  viscous  incompressible  fluid 
in  the  entrance  region  of  a  conduit  and  can  be  applied  to  an  annulus.   The 
method  of  linearization  was  first  originated  by  Targ  (1)  in  obtaining  the 
developing  velocity  profile  in  a  flat  duct  and  a  tube.   It  was  assumed  that 
the  fluid  entered  the  duct  with  a  uniform  velocity  and  that  the  transverse 
velocity  was  of  negligible  magnitude.   It  was  further  assumed  that  the  non- 
linear convective  term  could  be  replaced  by  the  product  of  the  initial  velo- 
city and  the  velocity  gradient  in  the  axial  dirvcttcn.   As  a  result,  the 


inertial  terms  of  the  Navier-Stokes  equation  were  linearized  and  a  linear 
momentum,  equation  obtained.   A  number  of  investigators  have  applied  this 
technique  of  Targ's  to  the  developing  flow  in  the  hydrodynamic  entrance 
region  of  annular  ducts;  among  them  are  Chang  and  Atabek  (2),  Roy  (3)  and 
Sparrow  and  Lin  (4).   The  latter  workers  solved  the  linearized  momentum 
equation  analytically  and  expressed  their  results  in  terms  of  modified  - 
Bassel's  functions. 

Another  important  approximate  method  of  solution  is  due  to  Langhaar 
(5),  who  obtained  steady  flow  patterns  in  the  transition  length  of  a  straight 
tube  by  means  of  a  linearizing  procedure.   He  introduced  the  assumption 
that  the  convective  terms  in  the  boundary  layer  equation  were  equal  to  the 
product  of  the  kinematic  viscosity  y,  a  parameter  $  which  is  a  function  only 
of  the  axial  distance,  and  the  axial  velocity  v  (r,z)  that  changed  with  the 
axial  and  radial  coordinates  z  and  r.   This  assumption  reduced  the  momentum 
equation  to  a  linear  form  and  thus  permitted  its  solution  by  analytical 
methods.   By  application  of  this  technique,  Sugino  (6)  and  Reynolds  et  al. 
(7)  have  obtained  analytical  velocity  profiles  for  the  developing  annular 
flow  problem.   A  third  method  based  on  a  series  solution,  has  been  used  by 
Murakawa  (8)  to  obtain  the  velocity  distribution,  the  pressure  drop,  and 
the  hydrodynamic  entry  length  in  an  annulus. 

Manohar  (9),  who  had  noted  that  considerable  difference  existed  between 
results  obtained  by  application  of  the  above  approximate  methods,  undertook 
an  exact  analysis  of  the  problem.   The  Navier-Stokes  equations  were  first 
simplified  under  the  usual  boundary  layer  assumptions  and  the  resulting  non- 
linear differential  equations  were  solved  by  an  implicit  finite-difference 
method.   The  idea  is  much  the  same  as  this  work  though  Manohar  did  not 


consider  the  heat  transfer  problem.   Comparisons  between  Manohar's  results 
and  those  of  this  work  have  been  made  and  will  be  discussed  in  later  sections. 

Very  little  experimental  work  has  been  done  on  the  developing  laminar 
flow  in  an  annulus.   Astill  (4)  has  provided  some  experimental  results  per- 
tinent to  the  analysis  of  Sparrow  and  Lin.   Velocity  profiles  were  measured 
using  air  as  the  fluid  and  were  compared  with  Sparrow  and  Lin's  theoretical 
results.   Though  agreement  between  experiment  and  analysis  was  considered 
reasonably  good,  it  was  reported  that  the  experimental  profiles  were  more 
skewed  than  the  analytical  profiles.   These  deviations  noted  in  the  neighbor- 
hood of  the  walls  were  explained  as  being  due  to  inevitable  difficulties  in 
measuring  small  velocities  near  solid  surfaces. 

Much  work  has  been  done  on  laminar-flow  annular  heat  transfer  during 
the  past  decade.   Most  of  the  published  papers  have  been  concerned  with  the 
thermal  entrance  region  where  a  fully  developed  velocity  profile  is  assumed 
at  the  inlet.   For  the  case  of  simultaneously  developing  velocity  and  tem- 
perature profiles  in  the  entrance  region  of  an  annulus,  however,  only  a  few 
publications  exist. 

Lundberg  et  al .  (10)  have  presented  a  theoretical  as  well  as  experimental 
analysis  for  the  thermal  problem  in  hydrodynamically  developed  flow.   This 
included  solution  of  an  eigenvalue  problem  with  four  different  sets  of  bound- 
ary conditions  listed  below: 

(i)   step  temperature  on  one  wall,  with  the  other  wall  maintained  at 
the  inlet  temperature, 

(ii)   step  heat  flux  on  one  wall  and  insulation  on  the  other, 

(iii)   step  temperature  on  one  wall,  with  the  other  v.all  remaining 
insulated, 


(iv)   step  heat  flux  on  one  wall  with  the  opposite  wall  maintained  at 
the  inlet  temperature. 

These  were  the  four  simplest  and  most  fundamental  conditions  that  could 
be  imposed  on  an  annulus.   The  linearity  and  homogeneity  of  the  energy  equa- 
tion permits  combination  of  a  variety  of  these  conditions.   This  technique 
of  superposition  can  be  best  illustrated  by  considering  some  of  the  rather 
complex  situations.   The  case  of  constant  heat  flux  on  one  wall  and  step 
temperature  on  the  other,  for  example,  can  be  considered  as  a  superposition 
of  conditions  (iii)  and  (iv)  while  the  case  of  constant  heat  fluxes  specified 
on  both  walls  required  superposition  of  condition  (ii)  twice,  one  imposed  on 
one  wall  and  the  second  on  the  other.   Viskanta  (11)  and  Hatton  and  Quarmby 
(12)  have  also  presented  solutions  to  some  of  these  same  problems.   Hsu  and 
Huang  (13)  considered  the  problem  from  a  slightly  different  viewpoint.   In- 
stead of  prescribed  wall  temperature  and  heat  fluxes ,  thermal  radiation  from 
the  walls  was  considered.   Newton's  law  of  cooling  was  employed  and  solutions 
were  obtained  for  cases  of  unilateral  (outer  wall  insulated)  and  bilateral 
radiative  transfer  and  also  for  the  case  where  radiative  heat  transfer  at 
the  wall  and  uniform  internal  heat  generation  take  place  simultaneously. 
Recently,  Hong  (14)  has  presented  a  study  of  the  thermal  entrance  problem 
with  fully-developed  velocity  profile  for  non-Newtonian  fluids  represented 
by  the  power-law  model.   Analytical  solutions  were  presented  for  different 
values  of  the  power-law  model  indices  and  of  the  inner  to  outer  radius 
ratio. 

There  are  two  theoretical  as  well  as  experimental  works  published  for 
the  case  of  simultaneously  developing  velocity  and  temperature.   Murakawa 
(8)  has  presented  numerical  and  experimental  results  for  the  case  of  constant 


wall  temperature  at  the  inner  wall  and  insulation  at  the  outer  wall.   Reynolds 
et  al.  (7)  have  presented  the  results  of  a  four  year  study  of  annular  heat 
transfer  to  Newtonian  fluids.   Included  in  their  study  is  a  complete  analysis 
of  the  annular  heat  transfer  problem  with  simultaneous  development  of  velo- 
city and  temperature  profiles.   The  Langhaar's  method  of  linearizing  approxi- 
mation was  adopted  and  the  problem  solved  by  an  integral  method.   Results  are 
tabulated  for  several  inner  to  outer  tube  radius  ratios  and  Prandtl  numbers. 
Experimental  measurements  were  made  using  air  (Pr  =  0.7)  and  agreement  between 
theory  and  experiment  was  considered  excellent. 

Shohet  (15),  using  a  numerical  approach,  solved  a  very  general  annular 
entry  problem  that  took  into  account  the  effect  of  a  magnetic  field  on  an 
electrically  conducting  fluid.   The  governing  equations  are  basically  the 
same  as  those  used  in  this  work  except  appearing  in  the  momentum  and  energy 
equations  are  several  more  terms  that  account  for  the  applied  magnetic  field. 
In  addition  to  the  Prandtl  number  and  the  inner  to  outer  radius  ratio,  a 
third  parameter,  the  Eckert  number,  comes  into  consideration  in  this  case. 
Results  were  presented  graphically  for  several  combinations  of  these  para- 
meters. 
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CHAPTER  III 


DEVELOPMENT  OF  THE  MATHEMATICAL  MODEL 
AND  FINITE  DIFFERENCE  EQUATIONS 


An  analysis  will  be  made  of  the  problem  of  simultaneous  development  of 
velocity  and  temperature  profiles  in  the  entrance  region  of  an  annulus.   The 
basic  governing  equations  which  are  used  in  analyzing  the  problem  will  be 
presented.   General  forms  of  these  equations  will  first  be  introduced  and 
then  reduced  by  means  of  simplifying  assumptions. 

1.   Mathematical  Statement  of_  Problem 

The  standard  equations  for  solving  this  type  of  problem  will  be  used, 
namely  the  continuity  equation,  the  momentum  equation  and  the  energy  equa- 
tion.  In  vector  form,  they  are  given  by 

Continuity: 

(1) 

Momentum: 


(2) 


Energy: 


pCv  Dc"  =  "  (^'q)  ~  T  (Jt)v(V'v)  "  (t:W)  (3) 


where  the  energy  equation  is  written  in  terms  of  the  temperature  T  of  the 
fluid. 

Cylindrical  polar  coordinates  r,  0  and  z  will  be  used,  with  the  z-axis 


lying  along  the  axis  of  the  tube  and  the  origin  of  the  coordinate  system  at 
the  center  of  the  inlet  cross  section.   The  annulus ,  together  with  the  co- 
ordinates and  the  corresponding  notation,  is  shown  in  Figure  1.   Being  axi- 
ally  symmetrical,  the  flow  will  be  independent  of  6.   Subject  to  the  follow- 
ing assumptions: 

(i)   The  fluid  is  incompressible  with  constant  physical  properties, 

(ii)   The  flow  is  steady  and  laminar, 

(iii)   There  is  no  viscous  dissipation, 

(iv)   Heat  conduction  in  the  z-direction  is  negligible, 

these  equations  reduce  to: 
Continuity: 


3v_ 
r  3r 

Momentum: 

z-component 

3v 


z     3p  ,   , 


Energy : 


(A) 


(    -)  + J]  (5) 


K  r3r    K  z  3z      3z  "  ^Lr  3rv  3r 
r-component 


3z 


"cplVrl?+'«lfl-  «ifc<«  £»  (7) 


By   utilizing  an   order-of-magnitude   analysis,    ' 


E 


o 

o 

c 


o 

8 
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the  Navier-Stokes   equations    considerably.      Following  his    arguments,    it    is 
possible   to   drop   the   term   for  the   molecular   transport   in   the   z    direction 
which   appears   in   the   z-component   of   the  momentum  equation.      It   is   also 
possible   to   drop   the   r-coraponent   of    cb.e   equation   completely   and   assume   that 
the   pressure  varies    only   in   the   z  direction.      As    a  result,    the  number   of 
unknowns   has   been   reduced  by   one.      Applying   these   assumptions,    the   governing 
equations   reduce   to: 

-  |-(rv  )   +  ~-  -  0  (8) 

r  3r       r         3z  . 

3v  8v  ,  ,      _      3v 

Examination  of  these  equations  shows  that  the  temperature  variable  is  pre- 
sent only  in  the  energy  equation.   Equations  (8)  and  (9)  can  then  be  solved 
simultaneously  to  give  the  velocity  profile.   Once  it  is  known,  the  tempera- 
ture profile  can  be  calculated  from  Equation  (10). 

2.   The_  Developing  Velocity  Profile 

Assuming  that  the  fluid  enters  the  annulus  with  a  uniform  velccitv  v 

o 

and  utilizing  the  fact  that  the  fluid  sticks  to  the  solid  surfaces  the  bound- 
ary conditions  for  Equations  (8)  and  (9)  may  be  written: 

z<0,       r,  <   x  <   r    .  v=v,       v  =  0 

—  i  —   —  o         z    o         r 

z  >  0,        r  =  ri'  vz  =  °'         Vr  =  °         (11) 

z>0,       r  ■  r  ,  v  «=  0,        v     •  0 


Examination  shows  that  there  are  two  equations  with  three  unknowns,  the  axial 
and  radial  velocities  and  the  pressure.   In  order  to  obtain  the  solution,  a 
third  equation  is  required.   As  is  customarily  done,  the  third  equation  is 
obtained  from  a  mass  balance  between  the  entrance  and  an  arbitrary  distance 
downstream.   The  resulting  equation  is, 

2tt  r  2    2 

/  /  °  pv  rdGdr  '=  pv  TT(r  -  r.)  (12) 

Z  O     O      1 

o  r. 

i 

This   equation   is   in   fact   an  integrated   form  of  the   equation   of   continuity. 
Discussion   of  the  use   of  this    apparently  non-independent  equation  may  be 
found   in   (16)    and    (17). 

For   convenience,    the    following  dimensionless   variables    are   introduced: 


V  =  v  /-y-  (13) 

rpr. 


r  2 


The  determining  equations  and  the  boundary  conditions  then  become 


(14) 
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r  /r.         ..   r  2 
/  °  X  URdR  «=  -^-[C-2)   -  1]  (16) 


Z  £  0, 

1    <   R   <   r    /r. , 
—      —    o     i 

U   =   1, 

V  =  0 

Z   >   0, 

R  =   1, 

U  =   0, 

V  =   0 

Z  >   0, 

R  =   r   /r. , 

O        l' 

U  =   0, 

V  =   0. 

(17) 


Note  that  the  only  parameter  that  appears  in  these  equations  and  the  boundary 

conditions  is  the  radius  ratio  r  /v..      Had  the  boundary  layer  assumptions 

not  been  introduced,  in  addition  to  the  ratio  r  /r. ,  the  Reynolds  number 
'  o  1*  J 

would  also  have  appeared  as  a  parameter. 

Thus,  while  considerable  mathematical  simplification  has  been  achieved, 
the  non-linear  character  of  the  Navier-Stokes  equation  has  been  preserved. 
There  remains  a  system  of  three  simultaneous  equations  to  be  solved  for  the 
three  unknowns  U,  V  and  P. 

A.   The  Finite-Difference  Technique 

Only  a  small  fraction  of  the  partial  differential  equations  generated 
in  engineering  problems  can  be  solved  by  formal  analytical  methods.   With 
the  non-linearity  of  the  momentum  equation  retained,  only  a  numerical  solu- 
tion can  be  provided.   The  finite-difference  technique  will  be  employed  in 
this  work.   By  using  difference  equations  in  place  of  the  derivatives  in  the 
partial  differential  equations,  the  problem  can  be  reduced  to  solving  sets 
of  algebraic  equations. 

To  elucidate  the  nature  of  this  method,  consider  a  mesh  network  imposed 
on  the  annulufl  as  shown  in  Figure  2.   The  mesh  covers  the  entire  cross  sec- 
tion of  the  annulus  and  continues  as  far  from  the  entrance  as  is  necessary* 
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Since  the  equation  is  independent  of  8,  it  is  understood  that  the  mesh  net- 
work applies  to  all  planes  normal  to  the  walls.   As  a  general  rule,  the  more 
the  mesh  is  refined,  the  smaller  the  truncation  error,  that  is,  the  smaller 
the  error  associated  with  determining  the  expressions  for  the  derivatives 
from  only  a  few  terms  of  the  Taylor  series. 

Equations  (14)  and  (15)  can  be  rewritten  as 


(18) 


V8R  +  U3Z     dZ    aR2    R  3R 

The  following  difference  quotients  will  be  used  for  the  derivatives  in  these 
equations : 

Continuity: 

3U  m   U(j+l,k+l)  -  U(j,k+1)  +  U(j+l,k)  -  U(j,k) 
az  2(AZ) 


(20) 


R(k+l)V(j+l,k+l)  -  R(k)V(j+l,k) 
AR 


3U  m   U(j-H,k)  -  U(j,k) 
3Z  AZ 

3U  m   U(j+l,k+l)  -  U(,j+l,k-l) 
3R  2(AR) 

32U  =  U(j+l,k+l)  -  2U(j+l,k)  +  U(j+l,k-l) 

2  2 

3R  (AR) 

d_P  _  P(j+1)  -  P(j) 
dZ"         AZ 


(21) 


The  difference  equations  for  Equations  (18)  and  (19)  at  a  typical  mesh  point 
(j ,k)  are  then 


R(k+l)V(j+l,k+l)  -  R(k)V(.j+l,k) 

AR 

+  YTaIO  ^R(k+1)tu(J+1>k+1)  "  U(j,k+1)]  +  R(k)[U(j+l,k)  -  U(j,k)]}  =  0 

(22) 

U(,j+l,k+l)-U(j+l)k-ll  U(j+l,k)-U(j,k) 

VU,K;         2(AR)         +  LU>k;       AZ 

_  P(j)-P(j+D  ,  U(,i+l,k+l)-2U(,i+l,k)+U(,j+l,k-l) 
*      iZ  (AR)2 

1  u(j+i,k+i)-u(.i+i,k-i) 

R(k)         2(AR)  K    ^ 


After  rearrangement  they  become 


v(j+1>k+1)  ■  R?Siy  v<j+1>k> 


-rr-  {R(k+l)[U(j+l,k+l)-U(j,k+l)] 


2(AZ)R(k+l) 
+  R(k)    [U(j+l,k)-U(j,k)]}  (24) 

C(j,k)U(j+l,k-l)   +  A(j,k)U(j+l,k)   +  3(j  ,k)U(j+l,k+l) 

+  P(,1+l)(AR)2   =   E(.>k)  (25) 
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2R(k) 

2 
A(j,k)  -  2  +"^f|-  U(j,k) 

(26) 
Ba,k)=fv(j,k)-1-^ 

Ea,k)«-^p(.j)+-^[u(j,k)]2 

Equations  (24)  and  (25)  can  be  applied  to  every  interior  point  k  for 
each  column  j.   Using  Equation  (25)  as  an  example,  the  following  equations 
can  be  written  for  each  k,  i.e.,  k  =  2,  3,  ,  n: 

C(j,2)U(j+l,l)+A(j,2)U(j+l,2)+B(j,2)U(j+l,3) 
+  P(j+l)(AR)2,E(j>2) 

C(j,3)U(j+l,2)+A(j,3)U(j+l,3)+B(j,3)U(j+l,4) 


»Ptf*l)W  ,K(3,3) 


(27) 


C(j,n)U(j+l,n-l)+A(j,n)U(j+l,n)+B(j,n)U(j+l,n+l) 

+  LU±IKM>!=E(J,n) 

In   these   equations,    U(j+l,l)    and   U(j+l,n+l)    represent    the   dlmetislo 
axial  velocity   at    the    inner  and   outer  wall   respectively    and  both   vanish   by 
the   second   and   third  boundary   conditions    of   Equation    (17).      The    same    La    true 
for  V(j+l,l)    and  V(j+l,n+l)    that   appear   in   Equation    (24). 
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The  material  balance  equation,  Equation  (16),  can  be  approximated  by 
the  trapezoidal  rule  for  column  j+1  as  follows: 


/r°/riURdR  =  AR  CU(J+l.DRq)+U(.i+l,2)R(2) 


U(j+1,2)R(2)+U(j+1,3)R(3)  , 
2 


U(j+l,n-l)R(n-l)+U(j+l,n)R(n) 
2  ' 


n-1 
AR  Z   U(j+l,k)R(k) 
k=2 


(28) 


where,  as  before,  both  U(j+l,l)  and  U(j+l,n)  vanish. 

The  finite  difference  approximations,  Equations  (20)  and  (21),  have 
been  employed  successfully  by  Shohet  (15)  in  obtaining  the  velocity  profile 
for  laminar  magnetohydrodynamic  flow  in  the  entrance  region  of  an  annulus. 
It  has  also  been  shown  (18) (19)  that  such  substitutions  produce  stable  dif- 
ference equations  for  any  selection  of  the  ratio  of  AR  and  AZ. 

B.   Solution  of  The  Velocity  Profile 

Equations  (27)  and  (28)  constitute  a  set  of  linear  algebraic  equations 
that  can  be  solved  simultaneously  to  give  the  pressure  and  the  n+1  values 
of  the  axial  velocity.   The  solution  was  obtained  by  a  trial  and  error  pro- 
cedure.  First,  a  reasonable  value  for  the  pressure  P(j+1)  in  Equation  (27) 
was  assumed.   The  resulting  set  of  n-1  equations  with  n-1  unknowns  were 
solved  by  the  Thomas  method  (20) .   Equation  (28)  was  then  used  as  an  equation 
of  constraint.   In  other  words,  it  was  used  to  determine  whether  the  correct 


values  of  U  were  obtained.   If  it  was  satisfied,  by  substituting  these  values 
of  D  into  Equation  (24),  the  radial  velocities  could  be  calculated.   Had  it 
not  been  satisfied,  another  value  for  the  pressure  would  have  been  assumed 
and  the  calculation  of  the  U's  repeated.   Examination  of  Equation  (27)  showed 
that  a  large  assumed  value  of  P(j+1)  gave  small  U's.   Thus,  if  the  left  hand 
side  of  Equation  (28)  was  greater  than  the  right  hand  side,  a  smaller  P(j+1) 
was  assumed  whereas  in  the  case  where  the  left  hand  side  was  smaller  than 
the  right  hand  side,  a  larger  P(j+1)  would  be  desired. 

To  obtain  the  velocity  profile  as  a  function  of  the  axial  and  radial 
distances,  calculations  were  started  at  the  j=l  column,  corresponding  to  the 
entrance  of  the  annulus.   Equations  similar  to  (25)  were  written  for  each 

interior  point  k=2,3,  ,n  while  Equation  (28)  was  applied  to  the  entire 

column.   The  resulting  system  of  n+1  algebraic  equations  were  then  solved 
according  to  the  above  procedure.   With  the  newly  found  axial  and  radial 
velocities  and  the  pressure  at  j=l,  calculations  could  proceed  to  the  j=2 
column.   In  this  manner,  integration  of  the  momentum  equation  was  completed 
from  j  to  j+1  and  it  was  possible  then  to  advance  column  by  column  along  the 
annulus  until  the  fully  developed  velocity  profile  was  obtained. 

The  calculations  were  performed  on  a  digital  computer  and  the  results 

are  presented  in  Figures  3,  4,  5  and  6  for  two  different  values  of  the  radius 

ratio,  namely  r  /r.  =  2.0  and  10.0.   The  flow  was  considered  fullv  developed 
J      o   l  '  ' 

when  the  velocity  at  j  and  j+1  remained  the  same  to  four  decimal  places. 
To  avoid  propagation  of  round-otf  errors  in  calculating  the  radial 
velocity  from  the  known  axial  velocity,  the  so-called  "biased  downward  velo- 
city" proposed  by  Shohet  (15)  was  introduced.   Contrary  to  Equation  (24)  in 
whicli  the  radial  velocity  was  calculated  in  an  upward  direction,  i.e.,  from 
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A/„ 
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the  inner  wall  to  the  outer  wall  of  the  annulus,  the  V's  were  also  calculated 
in  a  reversed  direction.   Thus,  by  substituting  the  following  difference 
quotients 

_3U  m   R(k-l)[U(j+l>k-l)-U(j,k-l)]+R(K)[U(,i+l>k)-U(,j>k)] 
3Z  2(AZ) 

(29) 

-l(vV\    =   R(k)V(j+l,k)  -  R(k-l)V(j+l,k-l) 
ZRK      J  AR 

into  Equation  (18) ,  a  slightly  different  form  of  the  continuity  equation  is 
obtained: 

R(k)V(j+l,k)  -  R(k-l)V(j+l,k-l) 


.  R(k-l) [U( j+l,k-l)-U( j ,k-l) ]+R(k) [U( j+1 ,k)-U(j ,k) ]  _ 
2(AZ) 


After  rearrangement,  it  becomes 


v(J+i,k-i)  -  R(k^)'k) 


2AZR(k-l) 

+  R(k) [U(j+l,k)-U(j ,k) ] }  (30) 

It  is  of  interest  to  note  that  by  using  the  same  values  of  the  axial 
velocity  calculated  from  Equation  (27)  ,  appreciably  different  results  were 
obtained  for  the  radial  velocity  from  Equations  (24)  and  (30)  respectively. 
Round-off  errors  were  accumulating  in  opposite  directions,  i.e.,  towards  the 
outer  wall  for  equation  (24)  and  towards  the  inner  wall  for  Equation  (30). 
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If  the  axial  velocities  were  exact,  then  within  computational  error,  Equation 
(24)  would  give  accurate  V  velocities  near  the  inner  wall  whereas  Equation 
(30)  would  give  accurate  V's  near  the  outer  wall.   To  minimize  these  round- 
off errors,  the  values  of  the  radial  velocities  calculated  both  from  Equa- 
tions (24)  and  (30)  were  averaged  for  all  but  the  interior  points  k=2  and 
k=n.   At  these  two  points,  instead  of  taking  the  average,  the  values  calcu- 
lated respectively  from  Equations  (24)  and  (30)  were  retained. 

As  was  mentioned  earlier,  calculation  of  the  axial  velocity  followed  a 
trial  and  error  procedure  in  which  the  macroscopic  material  balance,  Equation 
(28),  was  used  as  a  constraint.   An  allowable  error  of  approximately  ±0.01 
per  cent  was  imposed  upon  this  equation.   To  start  the  calculation,  the  first 
problem  to  be  confronted  is  the  value  of  the  axial  velocity  at  points  k=l, 
j  =  l  and  k=n+l,  j=l.   At  an  infinitesimal  distance  preceding  these  points, 
the  axial  velocity  is  one,  whereas  at  a  corresponding  infinitesimal  distance 
within  the  annulus,  the  axial  velocity  is  zero  on  the  walls.   In  other  words, 
the  flow  is  discontinuous  at  these  points  as  the  fluid  enters  the  annulus. 
While  there  was  no  information  to  determine  the  values  of  the  axial  velocity 
at  these  points,  an  average  value  of  0.5  was  assigned. 

Another  problem  is  the  proper  sizes  of  the  mesh  network  to  be  imposed 
upon  the  annulus.   In  general,  the  error  committed  in  replacing  the  deriva- 
tives of  the  boundary  layer  equations  by  difference  quotients  approached  zero 
as  the  mesh  was  refined.   Comparisons  could  be  made  by  using  several  different 
mesh  networks  and  the  most  appropriate  one  chosen.   Different  mesh  sizes  were 
used  along  the  annulus.   For  convenience,  they  are  listed  in  Table  1.   Re- 
sults will  be  compared  with  available  data  in  a  later  section. 


Table  1.   Mesh  sizes  for  the  continuity  and  momentum  equations. 
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Radius  ratio 

Mesh  number 

Vri 

of  R 

Z 

AZ 

AR 

2.0 

'51 

0.0 

0.00001 

0.02 

0.0015 

0.00005 

0.02 

0.0108 

0.0001 

0.02 

0.09 

0.0005 

0,02 

fully- 

developed 

10.0 

151 

0.0 

0.00001 

0.06 

0.0003 

0.00005 

0.06 

0.0024 

0.0001 

0.06 

0.009 

0.0002 

0.06 

0.045 

0.0005 

0.06 

0.108 

0.001 

0.06 

0.3 

0.002 

0.06 

1.14 

0.005 

0.06 

3.0 

0.01 

0.06 

7.8 

0.02 

0.06 

15.6 

0.05 

0.06 

fully- 

developed 

C.   The  Fully  Developed  Velocity  Profile 

The  laminar  flow  fully  developed  velocity  profile  in  an  annulus  has 
been  given  by  a  number  of  investigators.   The  derivation  is  presented  here 
for  two  reasons.   First,  the  definitions  of  the  dimensionless  variables  in 
this  work  are  different  from  those  of  previous  investigators  and  secondly 
the  equation  can  be  used  as  a  check  on  the  accuracy  of  the  numerical  solu- 
tion. 

The  appropriate  differential  equation  of  motion  for  fully  developed 
laminar  flow  is 


dZ   R  dRv   dR' 

and  the  boundary  conditions  are: 
R  =  1,         U  -  C 


(31) 


U«,  =  0        ,  (32) 

dUo, 


where  U^.  denotes  the  fully  developed  axial  velocity,  a  function  only  of  R. 
Integrating  Equation  (31)  twice  with  respect  to  R  yields 


(33) 


in  which  C   and  C  are  constants  of  integration.   Application  of  the  first 
and  third  conditions  of  Equation  (32)  gives 


°- -£#<»'  -  i-»Lln«  (34) 


whereas  application  of  the  second  and  third  condition::  loads  to 


2  dZ  "max 
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(35) 


Obviously,  Equations  (34)  and  (35)  must  yield  the  same  axial  velocity  for 
all  values  of  R.  Thus,  equating  the  two  equations  will  give  the  position 
where  maximum  velocity  occurs, 


2        _    (ro/r.)      r   1 

Rmov     =     9      In      (y      7Z     \  (36) 


velocity  defined  by 


(v  )  =  v 

z   avg         o 


2tt   r 
S     I  °  v  d0dr 

/     /  °rd0dr 
o     r. 


Equation    (37)    can  be  written   in   terms    of  dimensionless   variables    as 


2tt  r  /x. 
f     I  °     xu 


2tt   r  /r. 
/      /    °      XF 
o     1 


rn/r,- 

/  URdR 


2LV  Vli' 


(33) 


Substitution  of   Equations    (36)    and    (38)    into  Equation    (34)    results   in   the 
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following  expression    for   the   pressure   gradient, 


,_  (r  /r.)     -  1         r     2 

~  -  8/  C     °    *   ,,  ,     -   fc?)     -  1}  (39) 


Equation    (34)    can   then  be  written   as 


(r   /r   )*1  r     2 


(40) 


The  result  obtained  in  Equation  (40)  can  be  compared  with  the  numerical 
solution  obtained  far  from  the  entrance.   The  comparisons  are  made  in  Table 
2  and  it  caii  be  seen  that  agreement  is  excellent. 

D.  The  Pressure  Drop 

The  pressure  at  any  cross  section  along  the  annulus  was  obtained  simul- 
taneously with  the  axial  velocity  by  solving  Equations  (27)  and  (28) .   The 
pressure  drop,  which  is  due  to  friction  at  the  walls  and  the  change  of  mo- 
mentum between  the  entrance  and  any  downstream  cross  section,  could  then  be 
calculated  using  a  constant  dimensionless  pressure  of  10.0  at  the  inlet. 
With  an  allowable  error  of  approximately  ±0.01  percent  the  pressure  could 
be  calculated  to  the  fourth  significant  figure.   Results  are  presented  graph- 
ically in  Figures  7  and  8.   It  will  be  seen  in  a  later  section  that  the 
pressure  drop  is  not  only  a  function  of  the  dimensionless  distance  but  also 
depends  on  the  radius  ratio  and  decreases  with  increase  in  the  value  of  the 
ratio. 

E.  The  Entrance  Length 

An  estimate  of  the  hydrodynamic  entrance  length  may  be  made  on  th- 
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Table  2.   Analytical  and  numerical  fully-developed  velocity  profiles  for 


r  /r.  =  2.0 


J*< 


Radial  positions 

Fully-  dt 

R 

analytic; 

1.0 

0.0 

1.1 

0.603 

-1.2 

1.039 

1.3 

1.325 

1.4 

1.477 

1.471 

(R   ) 
max 

1.5075 

1.5 

1.502 

1.6 

1.411 

1.7 

1.210 

1.8 

0.904 

1.9 

0.5 

2.0 

0.0 

1.0 

0.0 

1.5 

0.558 

2.0 

0.924 

3.0 

, 

1.353 

4.0 

1.538 

4.64 

(R   ) 

max 

1.5673 

5.0 

1.558 

6.0 

1.450 

7.0 

1.230 

8.0 

0.911 

9.0 

0.499 

10.0 

0.0 

0.0 

0.6033 

1.0396 

1.3267 

1.4778 

1.5079 

1.5034 
1.4121 
1.2108 
0.9053 
0.5004 
0.0 


0.0 

0.5569 

0.9247 

1.352 

1.5378 

1.5672 

1.561 

1.4495 

1.2299 

0.9106 

0.4992 

0.0 


0.15 


Fig.    7.     Pressure  drop  in   the  hydrodynsmic  entrance   region,        /r.   =  2.0. 


10.0- 


2.0 


Fig.  8.   Pressure  drop  in  the  hydrodynamic  entrance  region,    It.    -    10.0. 


of   the  approach   of   the  numerical  axial  velocity   to  its    fully   developed  value. 


and   is    given  by 


,  r     -  In  R 

2[R2     -   1  -   {(-^)2   -    1}    =     .      m,a\] 
max  r ,  ln(r   /r. ) 


(41) 


ln(ro/r.)  \> 


(36) 


The  entrance  length,  Z  ,  is  then  defined  as  the  distance  at  which  the  maxi- 
mum velocity  attains  a  particular  percentage  of  the  maximum  fully  developed 
value,  that  is  , 


U(Z  ,  R   ) 
n    . =  desired  percentage.  (42) 


The  entrance  lengths  for  two  different  values  of  the  radius  ratio  are 
calculated  and  tabulated  in  Table  3  where  the  percentage  has  been  chosen  as 
99%. 

F.   Results  and  Discussion 

The  numerical  results  of  the  developing  velocity  profiles  are  shown  in 

Figures  3,  4,  5  and  6  corresponding  respectively  to  values  of  the  annulus 

r 
ratio  of  —  =  2.0  and  10.0.   The  good  agreement  between  the  analytical  and 

ri 
numerical  fully  developed  profiles  has  been  shown  previously  in  Table  2. 

To  further  demonstrate  the  adequacy  of  the  numerical  method  employed,  r> 
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of  this  work  will  be  compared  with  those  obtained  by  other  methods,  whenever 
they  are  available. 

Shown  in  Figures  9,  10  and  11  are  some  of  the  velocity  profiles  obtained 
by  a  number  of  investigators.   It  can  be  seen  that  for  both  values  of  the 
radius  ratio,  the  results  of  Man oh a r  (9),  who  likewise  employed  the  finite- 
difference  technique,  appear  to  deviate  slightly  from  those  of  this  work. 
His  results  are  considered  to  be  in  error  because  of  the  way  he  handled  the 
corner  conditions  at  the  point  of  entry.   The  usual  treatment  assumes  that 
the  fluid  approaches  the  entrance  with  a  uniform  velocity  profile,  including 
the  velocity  at  the  wall  positions.   Manohar,  on  the  other  hand,  assumed  a 
velocity  profile  where  the  velocities  were  taken  to  be  zero  on  the  boundaries 
and  unity  at  other  mesh  points.   In  other  words,  the  mass  flux  as  a  result 
of  his  assumption  of  the  inlet  conditions  is  less  than  usual  and  consequently 

leads  to  the  deviation  of  his  results. 

r 
The  calculated  results  by  Sugino  (6)  for  the  case  of  —  ■  2.0  at  R  =  1.5 

ri 
are  shown  in  Figure  9  for  comparison  and  his  velocity  profile  appears  to  de- 
velop much  faster.   The  results  obtained  by  Chang  and  Atabek  (2)  and  by 

r 
Sparrow  and  Lin  (4)  for  —  =  10.0  at  various  positions  are  compared  in  Figure 

,  ri 
11.   While  agreement  between  Sparrow  and  Lin's  results  and  that  of  this  work 

is  reasonably  good,  Chang  and  Atabek1 s  profiles  deviate  in  the  vicinity  of 
the  inlet. 

So  far  only  the  velocity  distributions  have  been  discussed.   The  changes 
in  these  profiles  along  the  annulus  have  an  important  effect  on  the  pressure 
drop.   As  can  be  seen  from  Figures  7  and  8,  the  pressure  drop  rises  sharply 
in  the  immediate  vicinity  of  the  inlet  and  grows  gradually  along  the  annulus 
until  it  finally  comes  to  a  constant  value  when  the   velocity  profile  becomes 
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fully-developed.   Figure  12  shows  the  pressure  drop  as  a  function  of  the 

axial  distance  and  the  radius  ratio,  where  the  axial  distance  is  taken  to 

be  Z'  =  for  convenience. 

r     2 

l 

A  comparison  between  Manohar's  results  (9)  and  those  of  this  work  is 

r 
provided  in  Figures  13  and  14  corresponding  to  —  =2.0  and  10.0  respectively. 

ri 
Again,  Manohar's  results  appear  to  be  too  small  and  such  deviation  is  be- 
lieved to  result  from  the  inaccuracy  of  his  velocity  profiles.   The  fully- 
developed  pressure  gradients  are  also  shown  in  these  Figures.   The  excellent 
agreement  between  analytical  and  numerical  results  once  again  substantiate 
the  accuracy  of  the  finite-difference  technique  employed.   In  Figure  15,  the 
pressure  drops  in  the  vicinity  of  the  inlet  are  plotted.   While  Sugino's 
results  (6)  seem  to  agree  with  that  of  this  work,  the  calculated  results  by 
Murakawa  (8)  appear  to  be  too  small  near  the  entrance. 

The  entrance  length  is  defined  by  Equation  (42)  in  a  previous  section. 
Most  of  the  published  work  provides  results  with  a  desired  percentage  of  99%. 
Since  the  methods  of  defining  the  dimensionless  axial  distance  by  various 
authors  are  different,  it  is  necessary  to  adopt  one  as  a  criterion  in  order 


to  compare  their  results.   As  before  the  dimensionless  quantity  that  takes 

z _ 

(Vri}  Uoph        <-* 


z              Z 
Into  account  the  radius  ratio,  Z'  =  « =  ,  is  employed. 


The  results  are  shown  in  Figure  16  with  Z'  plotted  against  the  recipro- 
r 
cal  of  the  radius  ratio  — .   Inspection  of  the  figure  reveals  that  the  Oft- 

ri 
trance  length  decreases  with  the  decreasing  radius  ratio.   It  also  shows  that 

the  results  obtained  by  Roy  (3)  are  almost  identical  to  those  by  Chang  and 

Atabek  (2)  whereas  the  results  of  Manohar  (9),  Reynolds  et  al .  (7)  and  this 
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Fig.  12.   Pressure  drop  in  the  hydrodynamic  entrance  region. 
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Fig.  13.  Comparison  of  pressure  drops  in  the  hydn 
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Fig.  14.   Comparison  of  pressure  drops  in  the  hydrodynamic  entrance  region, 
r 

/r-  =  10.0. 
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Fig.    15.      Pressure   drop    in    the    vicinity    of    the    Lnlet    for  —     =    2.0. 
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Fig.    16.      Comparison   of   entrance   lengths. 


work  are  approximately  the  same.   Sugino  (6)  concluded  in  his  work  that  the 
entrance  length  was  unaffected  by  the  radius  ratio  r  /r.  and  that  it  is 
approximately  0.02.   Based  on  the  calculations  of  this  and  other  work,  his 
conclusion  seems  in  error. 

3 .   The  Temperature  Profile 

The  problem  of  determining  the  developing  temperature  profiles  in  the 
entrance  region  of  an  annulus  will  be  presented  in  the  following  sections. 
With  a  numerical  method  of  solution,  any  kind  of  wall- temperature  variation 
or  heating  rate  can  be  handled  with  equal  ease.   Consideration  has  been 
restricted,  however,  to  two  distinct  limiting  cases. 

A.   Case  I:   constant  temperature  at  the  inner  wall,  outer  wall 
insulated. 

A.l.   Mathematical  statement  of  the  problem. 

The  coordinates  and  geometry  of  the  system  have  been  shown  in  Figure  1. 
The  fluid  flows  in  steady  laminar  motion  in  the  annulus  with  a  developing 
velocity  profile.   In  the  region  z  <  0,  the  fluid  and  both  walls  are  main- 
tained at  a  uniform  temperature  T  .   In  the  region  z  >  0,  a  constant  temper- 
ature, T  ?   T  ,  is  prescribed  at  the  inner  wall  and  the  outer  wall  is  in- 
sulated.  The  situation  is  shown  in  Figure  17a.   It  is  desired  to  find  the 
temperature  profile  and  the  variation  of  the  heat  transfer  coefficient  on 
the  inner  wall  with  axial  distance. 

Subject  to  the  limitations  mentioned  in  an  earlier  section  of  this 
chapter,  the  equation  of  energy  describing  the  problem  is 


■V*«H+,,.Si-MiiN&i  (10) 
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Fig.    17a.      Constant   temperature  at   the   inner  wall. 
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Fig.    17b.      Constant  heat   flux   at   the   inner  wall. 


with  boundary    conditions : 

z  ^  0,  r.    <   r   <   r    ,  T=T 

1  —       —     o  o 

z   >  0,  r  =  r    ,  T  =  1  (43) 

z    >  0,  r  =   r    ,  I1  =  0 

o  3r 

As  before,  introducing  Equation  (13)  and  the  following  definition  of 
the  dimension less  temperature, 


T  -  T 
T  -  = =?  (44) 


the  energy  equation  in  dimensionless  form  then  becomes 


3R    3Z   Pr  R  3RV  3Ry 


(45) 


where  Pr  is  the  Prandtl  number.   The  corresponding  boundary  conditions  are 
as  follows: 

Z  <  0,       l<R<r/r.,       T=0 
—  ——oi 

Z>0,       R=l,  T  •  1  (46) 

Z  >  0,       I  -  ro/r.,  H  -  0 

The  axial  and  radial  velocity  distributions  being  known,  Equation  (45) 
can  be  solved  readily  to  give  the  temperature  profile.   Unlike  the  velocity 
profiles,  solution  of  Equation  (45)  involves  no  trial  and  error  and  the 
procedure  is  straight-forward.   In  addition  to  the  radius  ratio  r  /r  ,  the 
Prandtl  number,  Pr,  appears  in  the  energy  equation.   The  effect  of  this 
dimensionless  quantity  will  be  discussed  later. 
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C(j,k)  =  -  -|v(j,k)  + 


2  VJ'  '        2R(k)Pr   Pr 


_2.iML_u(jjk) 


(50) 


B(j,k)  =  ~V(j,k)  - 


2*vj»w    2PrR(k)   Pr 


D(j,k)  =  -^-U(j,k)T(j,k) 


with  boundary  conditions 


Z  _<  0, 

1   <   R  <   r   /r.  , 

—        —      O      1 

T  =  0 

Z   >  0, 

R  =   1, 

T(j+l,l)    =   1 

Z   >   0, 

R  =  r   /r. , 

O        1 

T(j+l,n)    =  T(j+l,n+l) 

(51) 


The  quantities  T(j+l,n+l)  and  T(j+l,n)  represent  the  temperatures  at  the 
outer  wall  and  at  the  point  adjacent  to  the  wall  respectively. 

Equations  similar  to  Equation  (49)  can  be  written  for  every  interior 
point  k  for  a  particular  column  j+1.   Thus  for  k=2,  Equation  (49)  becomes 

C(j,2)T(j+l,l)+A(j,2)T(j+l,2)+B(j,2)T(j+l,3)  =  D(j  ,2) 
Application  of  the  second  boundary  condition  of  Equation  (51)  then  gives 

A(j,2)f(j+lJ2)+B(j,2)T(j+l,3)  =  D(  j  ,2)-C(  j  ,2)  (52) 

For  k=3,4,  ,  n-1,  the  following  equations  can  be  written 
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A- 2,      Solution  of  the  equation. 

To  obtain  the  temperature  profiles  by  solving  the  energy  equation,  the 
finite  difference  technique  with  the  Thomas  method  used  to  solve  the  result- 
ing algebraic  equations  will  again  be  employed.   The  following  finite  dif- 
ference approximations  are  introduced: 


il  =  T(j+1,K)  ~  T(j,k) 


II  =  T(j+l,k+l)  -  T(j+l,k-D  (      . 

3R  2(AR)  K      } 


9  T  _  T(j+l,k+l)  -  2T(j+l,k)  +  T(j+l,k-l) 

2  2 

3R  (AR) 

Substitution  of  these  approximations  into  the  energy  equation  yields 

T(j+l,k+l)-T(J+l,k-l)     (   .T(1+l,k)-T(1,k) 

_1  1_  T(j+l,k+l)-T(j+l,k-l)   _1  T^j+l,k+l)-2T(j+l,k)+T(j+l,k-l) 

(48) 


Pr  R(k)         2(AR)  Pr 


After  rearrangement,  the  energy  equation  takes  the  form 

C(j,k)T(j+l,k-l)+A(j,k)f(j+l,k)+3(j,k)T(j+l,k-H)=D(j,k)        (49) 
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C(j,3)T(j+l,2)+A(j,3)f(j+l,3)HB(jJ3)T(j+lJ4)  =  D(j,3) 

C(.i.,n-l)T(j+l,n-2)+A(j,n-l)f(j+l,n-l)+B(j,n-l)T(j+l,n)  -  D(j,n-1) 
For  k=n,  agaia,  some  rearrangement  is  necessary. 

C(j,n)T(j+l,n-l)+A(j,n)T(j+l,n)+B(j,n)T(j+l,n+l)=D(j,n) 
According  to  the  third  condition  of  Equation  (51) ,  it  can  be  rewritten  as 

C(j,n)T(j+l,n-l)  +  [A(j,n)+B(j,n)]T(j+l,n)=D(j,n)  (54) 

Equations  (52),  (53)  and  (54)  then  constitute  a  set  of  n-1  linear  alge- 
braic equations  with  the  same  number  of  unknowns.   They  can  be  solved  by  the 
Thomas  method  without  difficulty.   The  procedure  is  much  the  same  as  in  solv- 
ing the  momentum  and  continuity  equations  for  the  velocity  profile  except  in 
this  case  it  is  simpler  because  no  trial  and  error  is  involved.   Sample 
results  are  shown  graphically  in  Figures  18  and  19. 

The  problem  of  discontinuity  at  the  points  k=l,  j=l  and  k=n+l,  j=l  does 
not  exist  in  solving  the  energy  equation.  Examination  of  Equations  (52) , 
(53)  and  (54)  shows  that  no  information  on  the  values  of  T(l,l)  and  f(l,n+l) 
is  needed  to  start  the  calculation.  The  Thomas  method  is  applied  at  each 
column  and  the  values  of  the  radial  and  axial  velocities  are  directly  em- 
ployed in  these  equations.  Calculations  can  proceed  in  this  manner  until 
the  temperature  profile  no  longer  changes  with  axial  distance. 

A.  3.   The  bulk  temperature  and  the  Nusselt  number. 

The  determination  of  the  variation  of  the  Nusselt  number  on  the  inner 
surface  along-  the  length  of  the  annulus  is  the  main  purpose  of  this  work. 
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1.5 


Fig.  18.   Developing  temperature  profiles  for  Case  I,   °/r   =  2.0  and  P   ■  0.01. 
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J/a   =  H 
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But  before  the  expression  for  the  Nusselt  number  can  be  derived,  the  mixing- 
cup  or  bulk  temperature,  T,  ,  must  first  be  determined.   By  definition, 

2tt    r 
/     /  °V  Trd9dr 

\  '  s  r  Z  (55> 

/     /   °v  rdGdr 
r.    z 
o        1 

Introduction  of  the  dimensionless  quantities  defined  in  Equations  (13)  and 
(44)  gives 


2tt  r  /r. 
/   /  °   ^"UTRdedR 


b    2tt  r  /r. 

/   /  °   ^JRdSdR 


o   1 


r   2 

[(-V  -  i] 


r  /r. 
/  °  1UTRdR  ,  (56) 

1 


The  integral  in  Equation  (56)  can  be  evaluated  numerically  v/ith  the  known 
velocity  and  temperature  profiles.   Thus  for  column  j+1,  application  of  the 
trapezoidal  rule  leads  to 


b     r   2 


I     U(j+l,k)T(j+l,k)R(k)  (57) 


K-r)  -  i 


k=2 


in  which  U  vanishes  at  k=l  and  k=n+l. 
The  Nusselt  number  is  defined  by 


(58) 
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The  equivalent  diameter  D   and  the  heat  transfer  coefficient  are  defined  by 
the  following  equations : 


(59) 


(60) 


Therefore,  the  Nusselt  number  on  the  inner  surface  is 


-k(f*)       2r.(-*-l) 
v3r  r=r.     i  r. 


-   2(^  -1)  9R  R=1  (61) 


where  the  fact  that  the  dimensionless  inner  wall  temperature  is  unity  has 
been  utilized. 

To  evaluate  the  temperature  gradient  at  the  inner  wall,  the  following 
reasoning  is  adopted.   At  the  inner  wall,  the  dimensionless  equation  of 
energy,  Equation  (45)  ,  may  be  written  as 


3T  ,  .3  T 
"3R    ^ 
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(62) 


Substituting  Equation  (47)  into  the  above  equation,  the  following  difference 
equation  is  obtained  for  column  j+1. 


1_  T(,j+lt2)-T(.i+l,0)    T(j+l,2)-2T(j-KL,l)+T(j+l,0) 


R(l)       2AR 


(63) 


(AR) 


T(j+1,0)  is  the  temperature  of  a  hypothetical  point  inside  the  wall  as  shown 
in  the  sketch. 

T(J+1.2). 

T(j+1,D- 

T(j+1,0)- 


T  =  1 


inner  wall 


L 


Rearrangement  of  Equation  (63)  gives  T(j+1,0)  in  terms  of  T(j+1,2)  and 

f(j+l,l). 

Thus,  , 


T(j+1,0) 


(AR) 


2R(1)AR 


(AR) 


(AR) 


2    2R(1)AR 


2T(j+l,l)  -  (^|+1)f(j+l,2) 

'         AR 


(64) 


since  R(l)=l. 

The   temperature   gradient   at   the   inner  wall   is    then 
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(*L\  =  T(j+1,2)-T(j+1,0) 

V3R;R=1  2(AR) 


1      ra/.M    ox        2T(,j+l,l)-(   2+1)1(1+1,2). 
2(AR)ll"+1,Z;  AR  J 

2 


T(j+l,2)-T(,j+l,l) 


(65) 


Substitution  of  this    result  into   Equation   (61)    gives    finally 


2(_o  .  ^(j+i.D-Kj+i.a  -  (66) 

l  AR(1  -   T,)(l  -  -|) 


The  variation  of   the  Nusselt  number  with   the   axial  distance   is   plotted 
in  Figures    20   and   21.      Values   of  0.01,    1.0   and   10.0  have  been   considered    for 
the  Prandtl   number. 

B.      Case  II:      Constant  heat    flux   at   the  inner  wall  and  outer  wall 
insulated. 

B.l.      Statement    and  solution   of  the   problem. 

The   flow   conditions    for  the   constant   flux   case  are   the  same   as   in   the 
previous   section,   but    the  boundary   condition  on   the   inner  surface   is   differ- 
ent.     Instead  of  maintaining  the   inner  wall  at    constant   temperature,    a  uni- 
form heat    flux  is   prescribed  on   it    as   shown  in   Figure   17b.      The  boundary 
conditions   are   as    follows : 


z<0,  r.    <  r  <  r   ,  T=T 

—     '  i  —      —     o  o 

z   >  0,  r  =  r.,  -k~  =  q  (67) 

z   >  0,  r  =  r    ,  I1  =  0 
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A  convenient    form  for  the   dimensionless    temperature   in    this    case   is, 


T  -   T 

T  =  7~   ,  (68) 

qr./k 


where  q  is  the  constant  heat  flux  prescribed  at  the  inner  wall.   The  energy 
equation,  in  dimensionless  form,  then  becomes 


3Z   PrLR  Rv  dW 


and  the  corresponding  boundary  conditions  are 


(45) 


Z   <  0, 
Z   >   0, 

1    <   R  <   r   /r. , 
—        —      O       1 

R  =    1, 

T  =   0 
3R 

Z   >  0, 

R  =   r   /r., 

O        1 

11  _    0 
3R  "   °" 

(69) 


Introduction  of  the  difference  quotients  in  Equation  (47)  reduces  the 

energy  equation  to  equations  exactly  the  same  as  (49)  and  (50) .   At  the 

inner  wall,  combination  of  Equation  (65)  and  the  second  condition  of  (69) 
yields 


^T     m   T(j+l,2)-T(j+l,l) 
VW         R(l-A|) 


(70) 


and  at  the  outer  wall, 


_3T  m   T(j+l,n)-T(j+l,n+l) 
9R  AR 


T(j+l,n)  =  T(j+l,n+l)  (71) 

Therefore,  for  column  j+1,  the  following  equations  can  be  written  for  k=2, 
3,  ,  n: 


C(j,3)T(j+l,2)+A(j,3)T(j+l,3)+B(j,3)T(j+l,4)+D(j,3) 
(72) 


C(j,n)T(j+l,n-l)+[A(j,n)B(j,n)]T(j+l,n)=D(j,n) 

These  linear  algebraic  equations  have  been  solved  by  the  Thomas  method 

on  a  digital  computer  and  the  temperature  profiles  are  shown  in  Figures  22 

and  23  for  r  /r  =2.0  and  10.0  and  Pr=0.01,  1.0. 
o   1  ' 

B.2.   The  bulk  temperature  and  the  Nusselt  number 

Once  the  temperature  distribution  is  known,  the  bulk  temperature  can 
be  evaluated  numerically  from  Equation  (57) .   However,  for  the  case  of  uni- 
form heat  flux  prescribed  at  the  inner  wall,  a  simple  relation  exists  between 
the  bulk  temperature  and  the  axial  distance  from  the  inlet.   The  derivation 
is  as  follows: 

An  energy  balance  is  made  between  the  entrance  and  an  arbitrary  dis- 
tance z  downstream.   The  resulting  equation  is 
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T  -  T 
o 

qri 


Fiy.  22.   Developing  temperature  pi  CI,  °/r  =  2.0  and  Pr  =  0.01, 
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W, 
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2tt  r 
(27Tr.z)q  =  /  /  °  C  (T-T  )v  rd9dr  (73) 


which  states  that  the  heat  entering  through  the  walls  is  the  same  as  the 
difference  between  the  heat  transported  through  the  cross  section  at  z=0 
and  z=z.   T  is  a  reference  temperature  and  in  this  case  it  is  the  entering 
fluid  temperature.   By  introducing  the  corresponding  dimensionless  variables, 
Equation  (73)  becomes 


2ttZ    (       .   o 

Pr       \ 

o  1 


7      *J*4 

^  =  /  °   X  UTRdR  (74) 


Substitution  of  Equation  (56)  into  the  above  equation  then  yields  the  de- 
sired relation.   Thus, 


_1  .It  [(V  -  i] 

Pr   2  VVJ    iJ 


Pr[(^)2  -  1] 
i 


(75) 


The  bulk  temperature  at  any  arbitrary  distance  from  the  inlet  can  then  be 
evaluated  readily  from  Equation  (75)  without  knowing  the  temperature  dis- 
tribution across  the  annular  space  at  that  particular  distance. 

The  Nusselt  number,  in  terms  of  dimensionless  variables  for  this 


66 


problem  becomes 


-(H)R  =  1 


(76) 


i     T  -  T, 

w    b 


where  T  is  the  inner  wall  temperature  and  is  equal  to  T(j+l,l). 

The  bulk  temperature  and  Nusselt  number  are  plotted  versus  the  axial 

r 
distance  for  —  =  2.0  and  10.0  and  Pr  =  0.01,  1  and  10  respectively  in 

ri 
Figures  24  and  25. 

It  is  of  interest  to  note  that  the  bulk  temperatures  for  the  uniform 
heat  flux  case  calculated  from  Equations  (75)  and  (57)  differ  by  as  much 
as  10%  at  regions  very  close  to  the  inlet. 

C.   Results  and  Discussion 

The  developing  temperature  profile  in  the  entrance  region  of  an  annulus 
has  been  obtained  for  two  boundary  conditions  on  the  inner  surfaces,  namely 
(i)  constant  temperature  and  (ii)  constant  heat  flux.   In  both  cases  the 
outer  surface  is  insulated.   Sample  results  are  shown  graphically  in  Figures 
18,  19,  22  and  23  since  it  is  not  practical  to  give  temperature  distribution 
as  a  function  of  axial  distance  for  all  values  of  the  radius  ratio  and  the 
Prandtl  number. 

The  mesh  sizes  for  the  energy  equation  are  tabulated  in  Tables  4  and  5. 
Due  to  the  fact  that  the  variation  of  the  velocity  in  the  vicinity  of  the 
inlet  is  large,  very  small  increments  in  the  axial  direction  were  used  in 
the  momentum  and  continuity  equations.   These  were  gradually  increased  along 
the  annulus  until  a  final  value  of  1/20  was  reached  when  the  velocity  profile 


~i — i i r 


_,     ,       , r 


_L_i__._ 


o-    i 
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Table  4.   Mesh  sizes  for  the  energy  equation,  Pr  =  0.01. 


.ius  ratio 

Mesh  number 

r  /r. 

of  R 

Z 

AZ 

AR 

O   1 

0.0 

2.0 

51 

0.0015 

fully- 
developed 

0.0 

0.00001 
0.00005 

0.02 
0.02 

10.0 

151 

0.0003 

0.0024 

0.009 

0.045 

0.108 

0.24 

fully- 
developed 

0.00001 

0.00005 

0.0001 

0.0002 

0.0005 

0.001 

0.002 

0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
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Table   5.      Mesh   sizes    for    the   energy  equation,   Pr  =  1.0   and   10.0. 


lius  ratio 

Mesh  number 

r  /r. 

of  R 

Z 

AZ 

AR 

o   X 

0.0 

2.0 

51 

0.0015 
0.0108 
0.0906 
0.4206 
1.0206 
10.0 

0.0 

0.00006 

0.0003 

0.0006 

0.003 

0.006 

0.03 

0.02 
0.02 
0.02 
0.02 
0.02 
0.02 

10.0 

151 

0.0003 
0.0024 
/  0.009 
0.045 
0.108 
0.3 
1.14 
3.0 
7.8 

15.6 

51.0 

92.0 
380.0 

0.00006 

0.0003 

0.0006 

0.0012 

0.003 

0.006 

0.012 

0.03 

0.06 

0.12 

0.3 

0.5 

1.0 

0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
0.06 
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became  fully  developed.   For  the  energy  equation,  however,  increments  as 
large  as  six  times  the  ones  used  for  the  momentum  and  continuity  equations 
were  employed.   Since  the  temperature  distribution  was  determined  from  the 
known  velocity  profile,  it  is  believed  that  a  larger  axial  increment  for 
the  energy  equation  will  not  lead  to  significant  error. 

The  surface  and  bulk  temperature  variation  for  various  radius  ratios 
and  Prandtl  numbers  are  shown  in  Figures  26-29.   It  can  be  seen  that  the 
surface  (for  the  constant  flux  case)  and  bulk  temperatures  increase  very 
sharply  for  small  Prandtl  numbers,  say  0.01,  and  very  slowly  for  large 
Prandtl  numbers  such  as  10.0.   For  the  case  of  constant  heat  flux  at  the 
inner  wall,  the  bulk  temperature  is  a  linear  function  of  axial  distance  as 
defined  by  Equation  (75).   The  large  difference  between  the  surface  and  bulk 
temperatures  for  a  radius  ratio  of  10.0  is  expected  as  the  thermal  boundary 
layer  extends  over  such  a  small  fraction  of  the  cross  section  in  the  entrance 
region. 

The  variation  of  the  Nusselt  numbers  with  axial  distance  are  presented 
in  Figures  20,  21,  24  and  25.   A  convenient  form  of  the  axial  distance  de- 
fined as  Z  -  Z/Pr  has  been  used.   Results  of  the  asymptotic  Nusselt  number 
are  also  presented  in  Table  6,  along  with  analytical  results  obtained  by 
Lundberg  et  al.  (10).   The  maximum  deviation  is  shown  to  be  less  than  5  per 
cent  with  the  largest  deviation  occurring  when  the  Prandtl  number  is  0.01. 
This  trend  was  expected  since  for  fluids  with  very  low  Prandtl  number,  the 
temperature  profile  is  established  much  faster  than  the  velocity  profile 
and  slight  errcr  in  the  velocity  distribution  might  lead  to  significant  error 
in  the  temperature  profile.   In  contrast,  for  large  Prandtl  numbers,  Che 
velocity  profile  develops  more  rapidly  and  more  accurate  temperatui 
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Table   6.      Comparison  of   analytical   and  numerical  asymptotic  values   of   the 
Nusselt  number. 


Case   I:      Constant   temperature   at    the   inner  wall  and   outer  wall  insulated. 

Radius   ratio               Prandtl  number  Asymptotic  Nu-j_ 

r   /r.                                     Pr  Analytical              Numerical  %  error 
o     1 

2.0                                 0.01  5.738                        5.73  0.14 

1.0-  5.738                        5.738  0.0 

10.0  5.738                        5.742  0.07 

10.0                                 0.01  11.56                        11.07  4.24 

1.0  11.56                         11.574  0.12 

10.0  11.56                         11.587  0.23 


Case   II:      Constant  heat   flux   at   the   inner  wall  and   outer  wall   insulated. 


2.0  0.01 

1.0 

10.0 

10.0  0.01 

1.0 

10.0 


6.181 

6.145 

0.58 

6.181 

6.085 

1.54 

6.181 

6.189 

0.13 

11.90 

11.511 

3.27 

11.90 

11.94 

0.34 

11.90 

11.984 

0.71 
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could  be  obtained.   Both  the  velocity  and  temperature  distribution  develop 
at  an  approximately  similar  rate  for  Prandtl  numbers  near  unity.   There  are 
fluctuations  in  the  values  of  the  Nusselt  number  near  the  entrance.   However, 
they  smooth  down  quickly  as  the  profile  develops. 

Several  attempts  have  been  made  to  improve  the  calculated  Nusselt  num- 
ber.  By  means  of  a  linear  interpolation,  the  velocity  profile  for  a  mesh 
with  size  double  the  original  one  could  be  obtained.   The  temperature  profile 
and  thus  the  Nusselt  number  were  then  calculated  using  the  new  velocity  pro- 
file.  A  better  approximation  (21)  to  the  temperature  gradient  at  the  inner 
wall  has  also  been  attempted.   The  gradient  can  be  written  as 


3T 


~[4T(R+AR)  -  3T(R)  -  T(R+2AR)] 


(77) 


3R   2(AR) 

where  T(R)  represents  the  temperature  at  the  wall  and  f(R+AR)  and  f(R+2AR) 
represent  the  temperatures  at  two  consecutive  points  out  from  the  wall.   The 
situation  is  shown  in  the  sketch. 


/> 


■v 


T(R+2   R) 
T(R+   R) 


T(R) 


.«•-  T(R-    R) 


inner  wall 


Equation    (77)    can  be  written   in   terms    of   the   subscript   notation  as 


H  =   2CWI4^(J+1,2)    "    3T(J  +  1'1)    "   T(j+1,3)] 


(78) 


The  Nusselt  number  for  the  case  of  constant  wall  temperature  then  becomes 
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Nui  =  (T:   "  D^|^T(j+l,2)  -  3T(j+l,l)  -  T(j+1,3)]  (79) 


where  the  wall  temperature  T(j+l,l)  is  unity. 

For  the  case  of  constant  heat  flux  at  the  inner  wall,  application  of 
the  second  boundary  condition  of  Equation  (69)  yields 


Cf^R-l35  T(AR)"[45(j+1'2)  "  ^W.D  "  T(j+1,3)]  =  -1 


T(j+l,l)  =  ^-[2AR  +  4T(j+l,2)  -  T(j+1,3)]  (80) 

The  Nusselt  number  is  then 


Nu  =  2{-f  -   1)  ± —  (81) 

ri      T(j+l,l)  -  T.  / 


Had  the  temperature  of  a  hypothetical  point  inside  the  wall  been  taken 
into  consideration  in  place  of  T(R+2AR) ,  a  different  form  of  the  temperature 
gradient  would  have  resulted,  i.e., 

3R  -  2AR  i-^J'-**'   JlvJ,i'u/   *\j--!-/j  (82) 

Substitution  of  Equation  (64)  into  the  above  equation  gives,  after  rearrange- 
ment , 

_3T  =  (1+AR)  [T(j+1,2)  -  TQ+1,1)]  (g3) 

Rd-^f) 

The  wall  temperature  and  the  Nusselt  numbers  can  be  calculated  with  some 


79 


manipulation. 
Case  I: 


Nu  -  2(—  -  1)  (1+AR>  [T(.i+1,D  -  T(,j+1,2)] 
1   'ri  R(l-Tj(l-% 


(84) 


T(j+l,l)  -T(j+1,2)  +   (1+AR;  (85) 


Nu.  =  2(—  -  1) (86) 

1     ri      T(j+l,l)  "  T, 


Table  7  shows  a  comparison  of  the  wall  temperatures  and  Nusselt  numbers 
obtained  by  these  methods.   It  can  be  seen  that  all  this  rearrangement  does 
not  seem  to  provide  improvement  at  all. 

The  Nusselt  numbers  obtained  are  compared  with  the  analysis  made  by 
Reynolds  et  al.  (7)  in  Figures  30-32.   The  general  agreement  can  be  con- 
sidered to  be  reasonably  good  and  indicates  that  their  approximate  technique 
is  quite  accurate. 


Table  7.   Comparison  of  surface  temperature  and  Nusselt  number  obtained  from 
three  distinct  difference  equations,  r  /r.  =  10.0  and  Pr  =  0.01. 


0.24 
0.24 
0.24 

0.228 
0.228 
0.228 

0.222 
0.222 
0.222 


Surface  temperature 
Case  II 

2.046490 
2.046620 
2.043195 

2.021716 
2.021843 
2.018422 

2.008661 
2.008792 
2.005367 


Nusselt  number  Nu.: 
Case  I  Case  II 


11.0901 
11.0535 
11.7819 

11.0455 
11.0092 
11.7255 

11.0619 
11.0252 
11.7257 


11.5263 
11.5254 
11.5507 

11.5302 
11.5293 
11.5546 

11.5372 
11.5362 
11.5616 


The  asymptotic  Nusselt  numbers  are  11.56  for  Case  I  and  11.90  for  Case  II. 
#  values  calculated  from  Equations  (79),  (80)  and  (81). 
+  values  calculated  from  Equations  (84) ,  (85)  and  (86) . 
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CHAPTER  IV 


SUMMARY 


The  problem  of  simultaneous  development  of  velocity  and  temperature 
profiles  in  the  entrance  region  of  an  annulus  has  been  solved  by  the  use  of 
a  numerical  procedure.   The  system  consists  of  highly  nonlinear  partial  , 
differential  equations  which  cannot  be  solved  otherwise.   The  work  is  appli- 
cable to  all  sorts  of  similar  problems  whereas  approximate  methods,  such  as 
Targ's  linearization  technique,  are  limited  to  only  certain  cases.   For 
example,  it  would  be  a  relatively  simple  task  to  allow  for  different  boundary 
conditions  or  to  include  temperature-dependence  of  density  and  viscosity. 

All  papers  pertinent  to  the  work  have  been  discussed  in  Chapter  II. 
They  were  divided  into  three  categories,  namely  (i)  the  velocity  profile, 
(ii)  the  temperature  profile  and  (iii)  the  simultaneously  developing  velocity 
and  temperature  profiles.   Most  of  the  investigators  had  employed  either 
Targ's  linearization  method  or  Langhaar's  approximation  technique  to  linearize 
the  momentum  equation  and  obtained  analytical  solutions  to  the  problem  in 
terms  of  modified  Bessel's  functions.   Unlike  the  others,  Manohar  (9)  solved 
the  problem  numerically  by  an  implicit  finite-difference  method.   His  re- 
sults, however,  were  in  error  because  of  the  erroneous  way  he  treated  the 
entrance  conditions. 

The  momentum  and  continuity  equations  were  first  solved  simultaneously 
to  give  the  developing  velocity  profile  and  the  pressure  gradient.   Applica- 
tion of  the  finite-difference  technique  reduced  the  problem  to  a  set  of  lin- 
ear algebraic  equations  which  were  solved  by  the  Thomas  method  through  a 
trial  and  error  procedure.   To  render  the  problem  soluble,  a  third  equation, 


the  mass  balance,  was  introduced  as  a  constraint  equation.   Results  obtained 
were  presented  graphically  in  Chapter  III  and  extensive  comparisons  made 
with  available  data.   The  fully-developed  velocity  profile  was  compared  with 
analytical  results  in  Table  2  and  agreement  was  excellent.   The  developing 
velocity  profiles  obtained  by  Manohar  (9)  were  in  error  for  both  values  of 
the  radius  ratio.   While  agreement  between  Sparrow  and  Lin's  results  (4)  and 
that  of  this  work  was  reasonably  good  for  r  /r.  =  10.0,  Chang  and  Atabek's 
profiles  (2)  deviated  in  the  vicinity  of  the  inlet.   The  pressure  drop  re- 
sults of  Manohar  were  likewise  erroneous  due  to  the  inaccuracy  of  his  velo- 
city profile.   On  the  other  hand,  Sugino's  results  in  the  vicinity  of  the 
entrance  for  r  /r.  ■  2.0  seemed  to  agree  with  this  work.   A  comparison  of 
the  entrance  lengths  in  Figure  16  showed  that  the  results  obtained  by  Roy 
(3)  and  Chang  and  Atabek  (2)  were  almost  identical  whereas  the  results  of 
Manohar  (9),  Reynolds  et  al.  (7)  and  this  work  were  approximately  the  same. 
Sugino  (6)  concluded  in  his  work  that  the  entrance  length  was  unaffected  by 

the  radius  ratio  r  /r..   Bas 
o   1 

his  conclusion  is  in  error. 

The  values  of  the  velocity  profiles  obtained  were  substituted  into  the 
energy  equation  to  solve  for  the  temperature  profile.   Two  distinct  cases 
with  various  Prandtl  numbers  and  radius  ratios  were  considered: 

(i)   Constant  temperature  at  the  inner  wall,  outer  wall  insulated, 
(ii)   Constant  heat  flux  at  the  inner  wall,  outer  wall  insul.it  ed« 
Results  of  the  profiles  as  well  as  the  bulk  temperature  and  Nueeelt  number 
are  shown  graphically  in  Chapter  III.   Comparisons  were  made  with  results 
obtained  by  Reynolds  et  al.  (7)  for  the  constant  heat  flux  case  for 
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the  approximate  technique  used  by  Reynolds  yielded  quite  satisfactory  heat 
transfer  rates. 
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NOMENCLATURE 

Symbols 

C  ,  C         Specific  heat 
v»   p 

C  ,  C_         Arbitrary  constants 

D  Equivalent  diameter  defined  as  2(r  -  r.) 

e  o    1 

g  Gravitational  acceleration 

h.  Heat  transfer  coefficient  at  the  inner  wall,  defined  by 

Equation  (60) 

k  Thermal  conductivity 

N  .  Nusselt  number  at  the  inner  wall,  defined  by  Equation  (58) 

p  Static  pressure  of  fluid  in  arbitrary  cross  section 

p  Static  pressure  of  fluid  at  pipe  entrance 

P  Dimensionless  static  pressure  defined  as   "~ 

MC 

Pr  Prandtl  number  defined  as  — r-^ 

k 

q  Heat  flux 

r  Radius  from  pipe  center 

r.  Inner  radius  of  annular  space 

r  Outer  radius  of  annular  space 

R  Dimensionless  radius  defined  as  — 


1 


Dimension less  radius  of  maximum  velocity 

Time  variable 

Temperature 

Temperature  at  the  entrance  to  the  annulus 

Constant  temperature  at  inner  wall  for  case  I 


T  Wall  temperature 

T  Bulk  temperature  defined  by  Equation  (55) 

T  -  Tn  T  -   Tn 

T  Dimensionless  temperature  defined  as 


Ti  -  Tc  <•  Vk 

T  -  T 

T  Dimensionless  wall  temperature  defined  as  tt 

w  r  q  r./k 

T,  Dimensionless  bulk  temperature 

v 

U  Dimensionless  axial  velocity  defined  as  — 

J  v 

o 

Uro  Dimensionless  fully-developed  axial  velocity 

(U^,)   .        Dimensionless  maximum  fully-developed  axial  velocity 

v  Velocity 

v  Radial  velocity 

r  J 

v  Axial  velocity 

z 

V  Constant  velocity  at  the  inlet  of  the  annulus 

v 

V  Dimensionless  radial  velocity  defined  as  — , 


Axial  coordinate 

2 

Vorip 

Dimension less  axial  distance  from  the  inlet  defined  as  z/ 


Dimensionless  axial  distance  defined  as  Z/( —  -  1) 

i 

Dimensionless  entrance  length 

Dimensionless  axial  distance  defined  as  Z/Pr 


Greek  symbols 

jj  Viscosity 

p  Density 

6  Polar  coordinate 

t  Shear  stress 
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APPENDIX 

A.   The  Thomas  Method. 

Various  methods,  either  iterative  or  non-iterative,  for  solving  systems 
of  linear  algebraic  equations  have  been  developed  in  the  past.   In  obtaining 
the  numerical  solution  of  the  present  problem,  the  partial  differential 
equations  had  been  reduced,  by  means  of  the  finite-difference  technique  to  a 
set  of  linear  algebraic  equations  and  the  Thomas  method  was  applied.   The 
procedure  is  as  follows: 

Consider  a  system  of  algebraic  equations  with  the  form 


r  =  2,  3,  ,  n-1 

,  +  a  x  =  d 

n  n-1    n  n  n 

The  unknowns  x, ,  x_ ,  ....,  x  are  eliminated  from  these  equations  by  letting 

1   2  n 

Wl  =  al 

w  =a  -  c  q  .  r=2,3,...,n 


and 


8i =  Vwi 

d  -  c  g   , 
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x     =   g     -qxr1  r  =   1,    2,    >   n-1 

r        6r       nr  r+1  '      '  ' 

While  W,    q  and   g  are  computed  in   the  order  of  increasing  r,   the  calculation 
of  the  x's   is   in   the  order- of  decreasing  r. 

An  advantage   of   the  Thomas  method  over  the  usual  matrix  inversion 
method  or  Gauss   elimination  method  is    that   less   storage,  space   and   computing 
time  are   required. 
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Pwead  r0/ri,  L'L 
AR,  ZKAX,  Pr's 


Assure  P 


Calculate  axial 

and 
radial  velocities 


f   Satisfy  nnss  balance  1]~ 


Yes 


Calculate  temperature  \ 
profile,  bulk  tenperature 
and  Nussclt  number    / 


C    z  -  zmx  i\ 


No 


No 


Z  o  Z+AZ 


Fig.  33.     Flow  chart  for  the  comput«r  program. 
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B.   Computer  program. 


C     THE  SIMULTANEOUS  DEVELOPMENT  OF  VELOCITY  AND  TEMPERATURE 
C     PROFILES  IN  THE  ENTRANCE  REGION  OF  AN  ANNULUS 

DIMENSION  A(200) ,B (200) , C(200) ,D(200) ,W(200) ,G(200) 

DIMENSION  U(200),UN(200),V(200)-,VN(200),R(200),VM(200),VA(200) 

DIMENSION  AA(200) ,BB(200) ,CC(200) ,DD(200) ,DDD(200) ,T(200) ,TL(200) 

DIMENSION  GG(200) ,GGG(200) ,WW(200) ,WWW(200) ,T1(200) ,T2 (200) 

DIMENSION  TOl (2 , 200) , T02 (2 , 200) , TNI (2 , 200)  ,TN2 (2 , 200) ,PPR(5) 

DIMENSION  UU(200),W(200) 
100  FORMAT  (6F12.6) 
200  FORMAT  (I5,4F20.8) 
300  FORMAT  (2F20.8) 

400  FORMAT  (1H1.///3X, 'PL' ,10X, 'DELR' ,9X, 'DELZ' ,9X, 'Z' }8X, 'ITER') 
500  FORMAT  (/4X, 'THE  AXIAL  VELOCITY  PROFILE  UN(K)1/) 
510  FORMAT  (/4X,'THE  RADIAL  VELOCITY  PROFILE  VA(K) '/) 
600  FORMAT  (/4X,  'ANU1'  ,6X, '  ANU2  '  ,14X,  'TBU1'  ,16X,  'TBU2  '  ) 
700  FORMAT  (2FI0.4,2F20. 8) 
800  FORMAT  (4F20.8) 
900  FORMAT  (4F12.6,I5) 

910  FORMAT  (1H1, ///5X,' DELR', 8X,' DELZ' , 12X, ' Z' , 7X, 'PR' ) 
920  FORMAT  (/4X, 'THE  TEMPERATURE  PROFILE  FOR  CASE  1  T1(K)7) 
922  FORMAT  (/4X, 'THE  TEMPERATURE  PROFILE  FOR  CASE  2  T2(K)'/) 
930  FORMAT  (2F10.3) 
940  FORMAT  (10F13.8) 
C     CALL  THE  RATIO  RO/RI  BETA 

C     ASSIGN  VALUES  TO  BETA , DELR, DELZ ,ZMAX,  AND  P  WHERE  P  IS  THE  INLET 
C     PRESSURE  AND  ZMAX  THE  MAXIMUM  AXIAL  DISTANCE 

READ  (1,100)  BETA, DELR, DELZ, DELY, ZMAX, P 

READ  (1,930)  (PPR(K),K=1,2) 

L=0 

Z=0.0 

AN=(BETA-1.)/DELR 

N=AN 

NN=N+1 

MM=N-1 

MA-N-2 
C     AT  THE  POINTS  OF  CONTINUITY,  SET 

U(l)=0.5 

U(NN)=0.5 

T(l)=0.5 

T(NN)=0.0 

TL(1)=0.0 

TL(NN)=0.0 

T01(l,l)=0.5 

T01(2,l)=0.5 

TO2(l,l)=0.0 

TO2(2,l)=0.0 
C     THE  BOUNDARY  CONDITIONS  ARE 

DO  1  J=l,2 

DO  1  K=2,NN 

T01(J,K)=0.0 

TO2(J,K)=0.0 


1  CONTINUE 
DO  2  K=2,N 
U(R>1.0 
T(K)=0.0 

2  TL(K)=0.0 

THE  MATERIAL  BALANCE 

XX- (BETA-BETA- 1 . ) / (2 . *DELR) 

CALCULATION  OF  RADIAL  POSITIONS 

DO  3  K=1,NN 

V(K)=0.0 

AK-K 

3  R(K)=1.+(AK-1.)*DELR 
DN=DELR*DELR/DELZ 

ASSUME  A  REASONABLE  VALUE  FOR  THE  PRESSURE 
PN=9. 993286 
15  CONTINUE 
ITER=1 
MX— 5 
L=L+1 
Z=Z+DELZ 

IF(Z.GT.ZMAX)  GO  TO  90 
IF(L-l)  13,13,44 

13  DO  14  K-l.NN 
UU(K)=U(K) 

14  VV(K)-V(K) 
44  CONTINUE 

NAME  THE  PRANDTL  NUMBER  PR 

PR=0.01 

CALCULATION  OF  THE  COEFFICIENTS  IN  THE  ORDER  OF  INCREASING  R 

FOR  THE  TEMPERATURE  PROFILES 

DO  40  K=2,N 

AA(K>2./PR+DN*U(K) 

BB (K)-DELR*V (K) /2 . -DELR/ (2 . *R(K) *PR)-1 . /PR 

CC(K)-DELR/(2.*R(K)*PR)-DELR*V(K)/2.-l./PR 
40  DD(K)-DN*U(K)*T(K) 

WW(2)=AA(2) 

GG(2)=(DD(2)-CC(2))/WW(2) 

DO  50  K-3,MM 

WW(K)-AA(K)-CC(K)*BB(K-1)/WW(K-1) 

GG(K)=(DD(K)-CC(K)*GG(K-1))/WW(K) 

IF(GG(K).LT. 0.00000001)  GG(K)=0.0 
50  CONTINUE 

WW(N)=AA(N)+BB(N)-CC(N)*EB(N-1)/WW(N-1) 

GG(N)=(DD(N)-CC(N)*GG(N-1))/WW(N) 

DO  88  K-2.N 
83  DDD(K)*=DN*U(K)*TL(K) 

WWW(2)-CC(2)+AA(2) 

GGG(2)-(DDD(2)-CC(2)*DBLR*(1.-0.5*DELR))/WW(2) 

DO  99  X--3,MM 

MWW(K)-AA(K)-CC(K)*BB(K-1)/WWW(K-1) 

GGG(K)-(DDD(K)-CC(K)*(X;G(K-1))/WWW(K) 

I F (GGG (K) . LT . 0 . 00000001 )  GCC (K) =0 . 0 
99  CONTINUE 

W\%IW(N)-AA(N)+BB(N)-CC(N)*BB(N-1)/W\;W(N-1) 

GGG(N)=(DDU(N)-CC(N)*CGG(N-l))/i; 
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C      SOLUTION  OF  THE  VELOCITY  PROFILE  BY  THE  THOMAS  METHOD 
PL=PN 
DO  4  K=2,N 
A(K)=2.+DN*U(K) 
B(K)=V(K)*DELR/2.-l.-DELR/(2.*R(K)) 

4  C (K) =DELR/ ( 2 . *R (K) ) -1 . - V (K) *DELR/ 2 . 
22  CONTINUE 

DO  5  K=2,N 

5  D (K) =DN*P+DN*U (K) *U (K) -DN*PL 
W(2)=A(2) 

DO  6  K=3,N 
W(K)=A(K)-C(K)*B(K-1)/W(K-1) 

6  G(K)=(D(K)-C(K)*G(K-1))/W(K) 

C     THE  AXIAL  VELOCITIES  AT  WALLS  ARE  ZERO 

UN(1)=0.0 

UN(NN)=0.0 
C     THE  VELOCITY  IS  CALCULATED  IN  THE  ORDER  OF  DECREASING  R 

UN(N)=G(N) 

DO  7  K=1,MA 

I=N-K 

7  UN(I)=G(I)-UN(I+1)*B(I)/W(I) 
SUMU=0 . 0 

DO  8  K=2,N 

8  Sm-IU=SUMU+R(K)  *UN  (K) 

C     THE  PRESSURE  IS  DETERMINED  IN  SUCH  A  MANNER  THAT  AN  ERROR  OF 
C      ±0.01%  IS  ALLOWED  FOR  THE  MATERIAL  BALANCE 
IF(ABS(XX-SUMU)-0.04)  55,55,77 
77  CONTINUE 

IF (MX)  70,70,80 
70  IF(XX-SUMU)  80,55,75 

75  PL=PL-0.001 

ITER=ITER+1 

GO  TO  22 
80  IF(XX-SUMU)  85,55,95 
85  PL=PL40,0001 

ITER=ITER+1 

MX=5. 

GO  TO  22 
95  PL-FL-0. 00001 

ITER=ITER+1 

GO  TO  22 
55  CONTINUE 
C      CALCULATION  OF  RADIAL  VELOCITY  FROM  EQUATION  OF  CONTINUITY 
C      RADIAL  VELOCITIES  AT  WALLS  ARE  ZERO 

VN(1)=0.0 

TO(SN)=0.0 

DO  18  K=1,MM 
18  VN(K+1)=R(K) *VN(K) /R(K+1)- (DELR/ (2. *R(K+1) *DELZ) ) *(R(K+1) *(UN(K+1) 
1-U(K+1))+R(K)*(UN(K)-U(K))) 
C      THE  BIASED  DOWNWARD  VELOCITY  DISTRIBUTION 

VM(1)=0.0 

?M(NN)=0.0 

DO  19  K=1,MM 
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I=NK-K 
19  VM(I)=R(I+1)*VM(I+1)/R(I)+(DELR/(2.*R(I)*DELZ))*(R(I)*(UN(I)=U(I)) 
1+R(I+1)*(UN(I+1)-U(I+1))) 

C      CALCULATION  OF  THE  AVERAGE  RADIAL  VELOCITIES 

VA(1)=0.0 

VA(NN)=0.0 

VA(2)=VN(2) 

VA(N)=VM(N) 

DO  21  K=3,MM 
21  VA(K)=(VN(K)+VM(K))/2. 
C     THE  TEMPERATURE  PROFILES  * 

C      CASE  1  —  CONSTANT  TEMPERATURE  AT  THE  INNER  WALL  AND  OUTER  WALL 
C  INSULATED 

C     THE  WALL  TEMPERATURE  IS  ALWAYS  UNITY 

Tl(l)=1.0 
C     THE  TEMPERATURE  DISTRIBUTION  IS  CALCULATED  IN  THE  ORDER  OF 
C     DECREASING  R 

T1(NN)=GG(N) 

T1(N)=GG(N) 

DO  60  K=1,MA 

I=N-K 
60  T1(I)=GG(I)-T1(I+1)*BB(I)/WW(I) 

SUMT1O.0 

DO  66  K=2,N 
66  SUMT1=SUMT1+UN(K)*T1(K)*R(K) 
C     EVALUATION  OF  THE  AVERAGE  TEMPERATURE  AND  NUSSELT  NUMBER 

TBU1=2 . *DELR*SUMT1/(BETA*BETA-1 . ) 

ANU1-2 . * (BETA-1 . ) * (Tl(l) -Tl (2) ) / (DELR*(1. -TBU1) *(1 . -DELR/2 . ) ) 
C  CASE   2   —    CONSTANT  FLUX  AT   THE   INNER  WALL  AND  OUTER  WALL   INSULATED 

T2(NN)=GGG(N) 

T2(N)=GGG(N) 

DO   101   K=1,MA 

I=N-K 
101   T2(I)=GGG(I)-T2(I+1)*BB(I)/WWW(I) 

T2(1)=T2(2)+DELR*(1.-0.5*DELR) 
C     EVALUATION  OF  BULK  TEMPERATURE  AND  NUSSELT  NUMBER 

TBU2=2 . *Z/ (PR* (BETA*BETA-1 . ) ) 

ANU2=2.*(BETA-1.)/(T2(1)-TBU2) 

IF(L-6)  122,205,205 
122  CONTINUE 

DO  120  K=1,NN 

0(K)-UN(K) 

V(K)=VA(K) 

T(R)=T1(K) 
120  TL(K)=T2(K) 

P-PL 

PN=PL-0.00014 

GO  TO  15 
205  CONTINUE 

WRITE  (2.A00) 

WRITE  (3,900)  PL,DELR,DELZ,Z, tTEF 

WRITE  (3,500) 


DO  121  K=1,NN 

WRITE  (3,200)  K,UN(K),VA(K),T1(K),T2(K) 
121  CONTINUE 

WRITE  (3,600) 

WRITE    (3,700)    ANU1,ANU2,TBU1,TBU2 

TEMPERATURE   PROFILES    FOR  PR-1.0   AND   10.0 

DM=DELR*DELR/DELY 

DO  999  L=l,2 

PR=PPR(L) 

DO  210  K=2;N 

AA(K) =2 . /PR+DM*UU (K)  „ 

BB(K)=DELR*W(K)/2.-DELR/(2.*R(K)*PR)-l./PR 

CC (K) =DELR/ ( 2 . *R (K) *FR) -DELR* W (K) / 2 . - 1 . /PR 
210  DD(K)=DM*UU(K)*T01(L,K) 

WW(2)=AA(2) 

GG(2)=(DD(2)-CC(2))/WW(2) 

DO  220  K=3,MM 

WW(K)=AA(K)-CC(K)*BB(K-1)/WW(K-1) 

GG(K)=(DD(K)-CC(K)*GG(K-1))/WW(K) 

IF(GG(K) .LT. 0.00000001)  GG(K)=0.0 
220  CONTINUE 

WW(N)=AA(N)+BB(N)-CC(N)*BB(N-1)/WW(N-1) 

GG(N)  =  (DD(N)-CC(N)*GG(N--1))/WW(N) 

DO  230  K=2,N 
230  DDD(K)=DM*UU(K)*T02(L,K) 

WWW(2)=CC(2)+AA(2) 

GGG(2)=(DDD(2)-CC(2)*DELR*(1.-0.5*DELR))/WWW(2) 

DO  240  K=3,MM  t 

WWW(K)=AA(K)-CC(K)*BB(K-1)/WWW(K-1) 

GGG(K)=(DDD(K)-CC(K)*GGG(K-1))/WWW(K) 

IF(GGG(K).LT. 0.00000001)  GGG(K)=0.0 
240  CONTINUE 

WWW(N)=AA(N)+BB(N)-CC(N)*BB(N-1)/WWW(N-1) 

GGG(N)=(DDD(N)-CC(N)*GGG(N-1))/WWW(N) 

TN1(L,1)=1.0 

TN1(L,NN)=GG(N) 

TN1(L,N)=GG(N) 

DO  250  K=1,MA 

I=N-K 
250  TN1(L,I)=GG(I)-TN1(L,I+1)*BB(I)/WW(I) 

SUMTN1=0 . 0 

DO  260  K=2,N 
260  SUMTN1=SUMTN1+UN (K) *R(K) *TN1(L ,K) 

TBU1=2 . *DELR*SUMTN1/ (BETA*BETA-1. ) 

ANU>2.*(BETA-l.)*(TNl(L,l)-TNl(L,2))/(DELR*(l.-TBUl)*(l.-DELR/2.)) 

TN2(L,NN)=GGG(N) 

TN2(L,N)=GGG(N) 

DO  270  K=1,MA 

I=N-K 
270  TN2(L,I)=GGG(I)-TN2(L,I+1)*BB(I)/WWW(I) 

TN2(L,1)=TN2(L,2)+DELR*(1.-0.5*DELR) 


TBU2=2 . *Z/ (PR*(BETA*BETA-1 . ) ) 

ANU2=2  .  *  (BETA--1 . )  /  (TN2  (L ,  1) -TBU2) 

WRITE    (3,910) 

WRITE    (3,100)    DELR,DELY,Z,PR 

WRITE    (3,920) 

DC   275  K=1,NN 
275  WRITE    (3,940)    K,TN1(L,K)  ,TN2(L,K) 

WRITE    (3,600) 

WRITE    (3,700)    ANU1,ANU2,TBU1,TBU2 

DO   310  K=1,NN 

T01(L,K)=TN1(L,K) 
310   T02(L,K)=TN2(L,K) 

IF(Z-0. 00024)  999,999,280 
280  DO  290  K=1,NN 
290  PUNCH  300,TN1(L,K),TN2(L,K) 
999  CONTINUE 

DO  320  K=1,NN 

D(K)=UN(K) 

V(K)=VA(K) 

T(K)=T1(K) 
320  TL(K)=T2(K) 

L=0 

P=PL 

PN=PL-0. 00014 

IF(Z-0. 00024)  15,15,33 
33  DO  125  K=1,NN 
125  PUNCH  800,UN(K),VA(K),T1(K) ,T2(K) 

GO  TO  15  t 

90  STOP 

END 


ANNULAR  HEAT  TRANSFER  WITH 

SIMULTANEOUSLY  DEVELOPING  VELOCITY  AND 

TEMPERATURE  PROFILES 


by 

KHANH-DAN  HA 
.S.,    Cheng  Kung  University,    1966 


AN  ABSTRACT   OF  A  MASTER'S   THESIS 


submitted   in  partial   fulfillment   of  the 


requirements    for   the  degree 


MASTER  OF   SCIENCE 


Department   of   Chemical  Engineering 

KANSAS   STATE   UNIVERSITY 
Manhattan,    Kansas 

1969 


ABSTRACT 

A  numerical  analysis  is  made  of  the  problem  of  simultaneously  developing 
velocity  and  temperature  profiles  of  a  Newtonian  fluid  in  the  entrance  region 
of  an  annulus.  The  flow  is  assumed  laminar  and  the  fluid  to  possess  constant 
physical  properties.   Two  distinct  cases  are  considered  in  this  study: 

(i)    Constant  temperature  at  the  inner  surface,  outer  surface 
insulated, 

(ii)   Constant  heat  flux  at  the  inner  surface,  outer  surface 
insulated. 
Results  are  presented  for  Prandtl  numbers  of  0.01,  1.0  and  10.0  and  for 
values  of  the  ratio  of  outer  to  inner  radii  of  2.0  and  10.0.   Comparisons, 
when  possible,  have  been  made  with  results  obtained  by  other  investigators. 


